
%----------------------------------------------------------------
% Options
%----------------------------------------------------------------

% @#define alternative_taylor_rule=0 


%----------------------------------------------------------------
% Define variables and parameters
%----------------------------------------------------------------

var pstar_1 pstar_2 pstar_3 pstar_4 pstar_5
	psi_1 psi_2 psi_3 psi_4 psi_5
	phi_1 phi_2 phi_3 phi_4 phi_5
    pi p_1 p_2 p_3 p_4 p_5
    s s_1 s_2 s_3 s_4 s_5
    mctilde y c n tfp
    markup_1 markup_2 markup_3 markup_4 markup_5
    i a nu i_ann pi_ann y_nat y_nat_mis;

varexo  eps_i
        eps_a
        ;   

parameters  betta siggma varphi eta theta_1 theta_2 theta_3 theta_4 theta_5 phi_pi phi_y  
            rho_a rho_i sig_eps_i sig_eps_a  logpibar  
            logpbar_1 logpbar_2 logpbar_3 logpbar_4 logpbar_5
            logybar logybar_sss logmctildebar 
            logpstarbar_1 logpstarbar_2 logpstarbar_3 logpstarbar_4 logpstarbar_5 
            logsbar_1 logsbar_2 logsbar_3 logsbar_4 logsbar_5 logsbar 
            logmarkupbar_1 logmarkupbar_2 logmarkupbar_3 logmarkupbar_4 logmarkupbar_5 
            logpsibar_1 logpsibar_2 logpsibar_3 logpsibar_4 logpsibar_5 
            logphibar_1 logphibar_2 logphibar_3 logphibar_4 logphibar_5 
            logabar lognbar logtfpbar logibar 
            posmp;


%----------------------------------------------------------------
% Parametrization
%----------------------------------------------------------------

load paramfile; 
set_param_value('betta',betta)   
set_param_value('siggma',siggma)  
set_param_value('varphi',varphi)  
set_param_value('eta',eta)     
set_param_value('theta_1',theta_1)  
set_param_value('theta_2',theta_2)  
set_param_value('theta_3',theta_3)  
set_param_value('theta_4',theta_4)  
set_param_value('theta_5',theta_5)  
set_param_value('phi_pi',phi_pi)  
set_param_value('phi_y',phi_y)   

set_param_value('rho_a',rho_a)   
set_param_value('rho_i',rho_i)   

set_param_value('sig_eps_i',sig_eps_i)  
set_param_value('sig_eps_a',sig_eps_a)  

set_param_value('logpibar',log(pibar))   
set_param_value('logpbar_1',log(pbar_1))  
set_param_value('logpbar_2',log(pbar_2)) 
set_param_value('logpbar_3',log(pbar_3)) 
set_param_value('logpbar_4',log(pbar_4)) 
set_param_value('logpbar_5',log(pbar_5)) 
set_param_value('logybar',log(ybar)) 
set_param_value('logybar_sss',log(ybar_sss)) 
set_param_value('logmctildebar',log(mctildebar)) 
set_param_value('logpstarbar_1',log(pstarbar_1)) 
set_param_value('logpstarbar_2',log(pstarbar_2)) 
set_param_value('logpstarbar_3',log(pstarbar_3)) 
set_param_value('logpstarbar_4',log(pstarbar_4)) 
set_param_value('logpstarbar_5',log(pstarbar_5)) 
set_param_value('logsbar_1',log(sbar_1)) 
set_param_value('logsbar_2',log(sbar_2)) 
set_param_value('logsbar_3',log(sbar_3)) 
set_param_value('logsbar_4',log(sbar_4)) 
set_param_value('logsbar_5',log(sbar_5)) 
set_param_value('logsbar',log(sbar)) 
set_param_value('logmarkupbar_1',log(markupbar_1)) 
set_param_value('logmarkupbar_2',log(markupbar_2)) 
set_param_value('logmarkupbar_3',log(markupbar_3)) 
set_param_value('logmarkupbar_4',log(markupbar_4))
set_param_value('logmarkupbar_5',log(markupbar_5))  
set_param_value('logpsibar_1',log(psibar_1)) 
set_param_value('logpsibar_2',log(psibar_2)) 
set_param_value('logpsibar_3',log(psibar_3)) 
set_param_value('logpsibar_4',log(psibar_4)) 
set_param_value('logpsibar_5',log(psibar_5)) 
set_param_value('logphibar_1',log(phibar_1)) 
set_param_value('logphibar_2',log(phibar_2)) 
set_param_value('logphibar_3',log(phibar_3)) 
set_param_value('logphibar_4',log(phibar_4)) 
set_param_value('logphibar_5',log(phibar_5)) 
set_param_value('logabar',log(abar)) 
set_param_value('lognbar',log(nbar)) 
set_param_value('logtfpbar',log(tfpbar))  
set_param_value('logibar',log(ibar))   

set_param_value('posmp',posmp)   


%----------------------------------------------------------------
% Equilibrium conditions
%----------------------------------------------------------------

model; 
    exp(pstar_1) = eta/(eta-1) * exp(psi_1) / exp(phi_1);
    exp(pstar_2) = eta/(eta-1) * exp(psi_2) / exp(phi_2);
    exp(pstar_3) = eta/(eta-1) * exp(psi_3) / exp(phi_3);
    exp(pstar_4) = eta/(eta-1) * exp(psi_4) / exp(phi_4);
    exp(pstar_5) = eta/(eta-1) * exp(psi_5) / exp(phi_5);

    exp(psi_1) = exp(mctilde) * exp(c) ^ (-siggma) * exp(y) + theta_1 * betta * exp(pi(+1)) ^ (eta) * exp(psi_1(+1));
    exp(psi_2) = exp(mctilde) * exp(c) ^ (-siggma) * exp(y) + theta_2 * betta * exp(pi(+1)) ^ (eta) * exp(psi_2(+1));
    exp(psi_3) = exp(mctilde) * exp(c) ^ (-siggma) * exp(y) + theta_3 * betta * exp(pi(+1)) ^ (eta) * exp(psi_3(+1));
    exp(psi_4) = exp(mctilde) * exp(c) ^ (-siggma) * exp(y) + theta_4 * betta * exp(pi(+1)) ^ (eta) * exp(psi_4(+1));
    exp(psi_5) = exp(mctilde) * exp(c) ^ (-siggma) * exp(y) + theta_5 * betta * exp(pi(+1)) ^ (eta) * exp(psi_5(+1));

    exp(phi_1) = exp(c) ^ (-siggma) * exp(y) + theta_1 * betta * exp(pi(+1)) ^ (eta-1) * exp(phi_1(+1));
    exp(phi_2) = exp(c) ^ (-siggma) * exp(y) + theta_2 * betta * exp(pi(+1)) ^ (eta-1) * exp(phi_2(+1));
    exp(phi_3) = exp(c) ^ (-siggma) * exp(y) + theta_3 * betta * exp(pi(+1)) ^ (eta-1) * exp(phi_3(+1));
    exp(phi_4) = exp(c) ^ (-siggma) * exp(y) + theta_4 * betta * exp(pi(+1)) ^ (eta-1) * exp(phi_4(+1));
    exp(phi_5) = exp(c) ^ (-siggma) * exp(y) + theta_5 * betta * exp(pi(+1)) ^ (eta-1) * exp(phi_5(+1));

    1 = 1/5 * exp(p_1) ^ (1-eta) + 1/5 * exp(p_2) ^ (1-eta) + 1/5 * exp(p_3) ^ (1-eta) + 1/5 * exp(p_4) ^ (1-eta) + 1/5 * exp(p_5) ^ (1-eta);

    exp(p_1) = ( (1-theta_1) * exp(pstar_1) ^ (1-eta) + theta_1 * exp(pi) ^ (eta-1) * exp(p_1(-1)) ^ (1-eta) ) ^ (1/(1-eta));
    exp(p_2) = ( (1-theta_2) * exp(pstar_2) ^ (1-eta) + theta_2 * exp(pi) ^ (eta-1) * exp(p_2(-1)) ^ (1-eta) ) ^ (1/(1-eta));
    exp(p_3) = ( (1-theta_3) * exp(pstar_3) ^ (1-eta) + theta_3 * exp(pi) ^ (eta-1) * exp(p_3(-1)) ^ (1-eta) ) ^ (1/(1-eta));
    exp(p_4) = ( (1-theta_4) * exp(pstar_4) ^ (1-eta) + theta_4 * exp(pi) ^ (eta-1) * exp(p_4(-1)) ^ (1-eta) ) ^ (1/(1-eta));
    exp(p_5) = ( (1-theta_5) * exp(pstar_5) ^ (1-eta) + theta_5 * exp(pi) ^ (eta-1) * exp(p_5(-1)) ^ (1-eta) ) ^ (1/(1-eta));

    exp(s) = 1/5 * exp(s_1) + 1/5 * exp(s_2) + 1/5 * exp(s_3) + 1/5 * exp(s_4) + 1/5 * exp(s_5);

    exp(s_1) = (1-theta_1) * exp(pstar_1) ^ (-eta) + theta_1 * exp(pi) ^ (eta) * exp(s_1(-1));
    exp(s_2) = (1-theta_2) * exp(pstar_2) ^ (-eta) + theta_2 * exp(pi) ^ (eta) * exp(s_2(-1));
    exp(s_3) = (1-theta_3) * exp(pstar_3) ^ (-eta) + theta_3 * exp(pi) ^ (eta) * exp(s_3(-1));
    exp(s_4) = (1-theta_4) * exp(pstar_4) ^ (-eta) + theta_4 * exp(pi) ^ (eta) * exp(s_4(-1));
    exp(s_5) = (1-theta_5) * exp(pstar_5) ^ (-eta) + theta_5 * exp(pi) ^ (eta) * exp(s_5(-1));

    exp(mctilde) =  exp(s) ^ varphi * exp(c) ^ siggma * exp(y) ^ varphi * exp(a) ^ (-(1+varphi));

    exp(c) ^ (-siggma) = betta * exp(c(+1)) ^ (-siggma) * exp(i) / exp(pi(+1));

    exp(y) = exp(c);

    exp(n) = exp(y) / exp(a) * exp(s);

    exp(tfp) = exp(a) * exp(s) ^ (-1);

    exp(markup_1) = exp(p_1) / exp(mctilde);
    exp(markup_2) = exp(p_2) / exp(mctilde);
    exp(markup_3) = exp(p_3) / exp(mctilde);
    exp(markup_4) = exp(p_4) / exp(mctilde);
    exp(markup_5) = exp(p_5) / exp(mctilde);
 

    @#if alternative_taylor_rule==0
        exp(i) = rho_i * exp(i(-1)) + (1 - rho_i) * (exp(logibar) * (exp(pi) / exp(logpibar)) ^ phi_pi * (exp(y) / exp(y_nat)) ^ phi_y) + nu;
    @#endif
    @#if alternative_taylor_rule==1
        exp(i) = rho_i * exp(i(-1)) + (1 - rho_i) * (exp(logibar) * (exp(pi) / exp(logpibar)) ^ phi_pi * (exp(y) / exp(y_nat_mis)) ^ phi_y) + nu;
    @#endif

    nu = posmp * sig_eps_i * eps_i;

    a  = rho_a * a(-1) + sig_eps_a * eps_a;

    i_ann = 4*i;
    pi_ann = 4*pi;

    exp(y_nat)^(-(siggma+varphi))*exp(a)^(1+varphi) = eta/(eta-1);
    exp(y_nat_mis)^(-(siggma+varphi))*exp(tfp)^(1+varphi) = eta/(eta-1);
end; 

 
%----------------------------------------------------------------
%  steady states
%---------------------------------------------------------------- 

initval;
p_1 = logpbar_1;
p_2 = logpbar_2;
p_3 = logpbar_3;
p_4 = logpbar_4;
p_5 = logpbar_5;
pi = logpibar;
y = logybar;
n = lognbar;
mctilde = logmctildebar;
pstar_1 = logpstarbar_1;
pstar_2 = logpstarbar_2;
pstar_3 = logpstarbar_3;
pstar_4 = logpstarbar_4;
pstar_5 = logpstarbar_5;
psi_1 = logpsibar_1;
psi_2 = logpsibar_2;
psi_3 = logpsibar_3;
psi_4 = logpsibar_4;
psi_5 = logpsibar_5;
phi_1 = logphibar_1;
phi_2 = logphibar_2;
phi_3 = logphibar_3;
phi_4 = logphibar_4;
phi_5 = logphibar_5;
s_1 = logsbar_1;
s_2 = logsbar_2;
s_3 = logsbar_3;
s_4 = logsbar_4;
s_5 = logsbar_5;
markup_1 = logmarkupbar_1;
markup_2 = logmarkupbar_2;
markup_3 = logmarkupbar_3;
markup_4 = logmarkupbar_4;
markup_5 = logmarkupbar_5;
c = logybar;
a = logabar;
i = logibar;
s = logsbar;
eps_i = 0;
tfp = logtfpbar;
eps_a = 0;
eps_i = 0;
nu = 0;
y_nat = logybar;
% i_ann = logi_annbar;
% pi_ann = pi_annbar;
end; 

steady;
check;


%----------------------------------------------------------------
% Shocks
%----------------------------------------------------------------

shocks;
    var eps_i = 1;   
    var eps_a = 1;   
end; 



%----------------------------------------------------------------
% Get variable positions
%----------------------------------------------------------------

% fun_var_pos


%----------------------------------------------------------------
% Solve
%----------------------------------------------------------------

if order==1
    stoch_simul(order=1, irf=17, irf_shocks=(eps_i, eps_a), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann n mctilde s s_1 s_2 s_5 pstar_1 pstar_5 p_1 p_5 markup_1 markup_5; 
elseif order==2
    stoch_simul(order=2, pruning, irf=17, irf_shocks=(eps_i), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann n mctilde s s_1 s_2 s_5 pstar_1 pstar_5 p_1 p_5 markup_1 markup_2 markup_3 markup_4 markup_5; 
elseif order==3 
    stoch_simul(order=3, irf=17, irf_shocks=(eps_i, eps_a), nomoments, nodecomposition, nocorr, nofunctions, noprint) pi tfp y i_ann markup_1 n mctilde; 
end
